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Abstract Dense granular clusters often behave like macro- 
particles. We address this interesting phenomenon in a model 
system of inelastically colliding hard disks inside a circular 
box, driven by a thermal wall at zero gravity. Molecular dy- 
namics simulations show a close -packed cluster of an almost 
circular shape, weakly fluctuating in space and isolated from 
the driving wall by a low-density gas. The density profile of 
the system agrees very well with the azimuthally symmetric 
solution of granular hydrostatic equations employing consti- 
tutive relations by Grossman et al, whereas the widely used 
Enskog-type constitutive relations show poor accuracy. We 
find that fluctuations of the center of mass of the system are 
Gaussian. This suggests an effective Langevin description 
in terms of a macro-particle, confined by a harmonic poten- 
tial and driven by a delta-correlated noise. Surprisingly, the 
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fluctuations persist when increasing the number of particles 
in the system. 
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1 Introduction 

Despite much progress in the last two decades, modeling of 
granular flow remains a challenge for physicists and engi- 
neers |6j; [15Q . The simplest case for a first-principle mod- 
eling seems to be a granular gas, or rapid granular flow: a 
flow dominated by binary particle collisions Q; [Till . Here 
one can use the model of inelastically colliding hard spheres 
to develop kinetic and hydrodynamic descriptions. The sim- 
plest version of this model involves binary collisions with a 
constant coefficient of normal restitution £: 
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where primed quantities stand for post-collisional velocities, 
and e, 7 = (r,- - r y ) / |r,- - r,-|. 

Assuming molecular chaos allows to use the Boltzmann 
or Enskog kinetic equation, properly generalized to account 
for the inelasticity of the particle collisions. Systematic deriva- 
tions of the hydrodynamic equations from the Boltzmann or 
Enskog equation Il3l; l2ll ; l2ql use an expansion in the Knudsen 
number and, in some versions of theory, in inelasticity [ 26] . 
A finite inelasticity immediately brings complications II ill . 
Correlations between particles, developing already at a mod- 
erate inelasticity, may invalidate the molecular chaos hy- 
pothesis I2l;l25l; l27ll . The normal stress difference and devi- 
ations of the particle velocity distribution from the Maxwell 
distribution also become important for inelastic collisions. 
As a result, the Navier-Stokes granular hydrodynamics should 
not be expected to be quantitatively accurate beyond the limit 
of small inelasticity. 
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A further set of difficulties for hydrodynamics arise at 
large densities. In mono-disperse systems strong correlations 
appear, already for elastic hard spheres, at the disorder-order 
transition. As multiple thermodynamic phases may coexist 
there [7], a general continuum theory of hard sphere fluids, 
which has yet to be developed, must include, in addition 
to the hydrodynamic fields, an order-parameter field. Fur- 
thermore, even on a specified branch of the thermodynamic 
phase diagram, we do not have first-principle constitutive 
relations (CRs): the equation of state and transport coeffi- 
cients. Last but not least, kinetic theory of a finite-density 
gas of elastic hard spheres in two dimensions has notorious 
difficulties related to the long-time tail in the velocity pair 
autocorrelation function 12811 . 

It has been argued recently that discrete particle noise 
may play a dramatic role in rapid granular flow lfH: [l2l ; l23ll . 
A new challenge for theory is a quantitative account of this 
noise. A promising approach at small or moderate densi- 
ties is "Fluctuating Granular Hydrodynamics": a Langevin- 
type theory that takes into account the discrete particle noise 
by adding delta-correlated noise terms in the momentum 
and energy equations, in the spirit of the Fluctuating Hy- 
drodynamics by Landau and Lifshitz II 1 8fl . There is a recent 
progress in this area in the case of small densities 0]. For 
high densities such a theory is presently beyond reach. 

This work addresses granular hydrodynamics and fluc- 
tuations in a simple two-dimensional granular system under 
conditions when existing hydrodynamic descriptions ll3l:l2ll 
[26h break down because of large density, not large inelas- 
ticity. In view of the difficulties mentioned above, attempts 
of a first-principle description of hydrodynamics and fluc- 
tuations should give way here to more practical, empiric or 
semi-empiric, approaches. One such approach to a hydrody- 
namic (or, rather, hydrostatic, as no mean flow is present) 
description was suggested in 1997 by Grossman et al. fl3ll . 
The present work puts it into a test in an extreme case when 
macro-particles (granular clusters with the maximum den- 
sity close to the hexagonal close packing) form. The model 
system we are dealing with was first introduced by Esipov 
and Poschel 19J]. It is an assembly of N ;§> 1 identical disks 
of mass in, diameter d and coefficient of normal restitution 
£, placed inside a circular box of radius R at zero gravity. 
The circular wall of the box is kept at constant temperature 
7b. We measure, using molecular dynamics (MD) simula- 
tions, the radial density profiles of the system, including the 
close-packed part. Furthermore, we solve numerically a set 
of granular hydrostatic equations which employ the CRs by 
Grossman et al. and show that, in a wide range of param- 
eters, there is good agreement between the two. We also 
show that, for the same setting, the Enskog-type CRs lfl4ll 
perform poorly. Finally, we investigate in some detail the 
fluctuations of the macro-particle's position, by measuring 
the radial probability distribution function (PDF) of the cen- 
ter of mass of the system. These fluctuations turn out to be 
Gaussian, which suggests an effective Langevin description 
of the system in terms of a macro-particle, confined by a har- 
monic potential and driven by a white noise. Surprisingly, 




Fig.l Snapshots of the system for N = 1716 (top left), 2761 (top right) 
and N = 5055 (bottom left). Also shown are magnifications of the in- 
dicated areas. 



the fluctuations persist as the number of particles in the sys- 
tem is increased. 



2 MD simulations 

We employed a standard event-driven algorithm and thermal 
wall implementation [24]. We put m = d = Tq = I and fixed 
R = 100 and e = 0.888, while the total number of particles ./V 
served as the control parameter in the simulations. We were 
mostly interested in a hydrodynamic (low Knudsen number) 
regime, when the mean free path is small compared to the 
system size. This requires N 3> R/d = 100. For N in the 
range of a few hundred, we observed a dilute granular gas 
with an increased density in the center of the box. The clus- 
tering in the center becomes more pronounced as N grows. 
The clustering can be easily explained in the hydrodynam- 
ics language: because of the inelastic collisions the granular 
temperature goes down as one moves away from the wall to- 
ward the center of the box. Combined with the constancy of 
the pressure throughout the system, this causes an increased 
particle density at the center. As N increases, the particle 
density in the center approaches the hexagonal close packing 
value n c = 2/(\^3d 2 ). Figure Q] shows snapshots of the sys- 
tem for three different, but sufficiently large, values of N. A 
perfect hexagonal packing is apparent. Movies of these sim- 
ulations show that the cluster position fluctuates around the 
center of the box, while the cluster shape fluctuates around 
a circular shape. 

Our diagnostics focused on two radial distributions: the 
number density of the particles n(r) and the PDF of a ra- 



3 



dial position of the center of mass of the system P(r m ), see 
below. As we are interested in steady-state distributions, we 
disregarded initial transients. 



3 Hydrostatic theory 

As the fluctuations are relatively weak, it is natural to start 
with a purely hydrodynamic description. In the absence of 
time-dependence and for a zero mean flow this is essentially 
a hydrostatic theory. It operates only with (time-independent) 
granular density «(r), temperature T(r) and pressure p(r). 
The energy input at the thermal wall is balanced by dissipa- 
tion due to inter-particle collisions, and one can employ the 
momentum and energy balance equations: 



p = const, V- (jcVr) =/. 



(2) 



Here K is the thermal conductivity and / is the rate of energy 
loss by collisions. Note that the heat flux, entering the ther- 
mal balance in Eq. (0, does not include an inelastic term, 
proportional to the density gradient IT3l: l2ll : l25l . In the nearly 
elastic limit 1 — £ <C 1, that we are interested in, this term can 
be neglected. The boundary condition at the thermal wall is 



T(r = R,$) = T , 



(3) 



where r and <j> are polar coordinates with the origin at the 
center of the box. 

To proceed from here, we need CRs: an equation of state 
p = p(n,T) and relations for K and / in terms of n and T. 
As we attempt to describe close-packed clusters, the stan- 
dard techniques, based on the Boltzmann or Enskog equa- 
tions, are inapplicable. Grossman et al. |[l3ll suggested a set 
of semi-empiric CRs in two dimensions, that are valid for 
all densities, all the way to hexagonal close packing. Their 
approach ignores possible phase coexistence beyond the dis- 
order-order transition and assumes that the whole system is 
on the thermodynamic branch extending to the hexagonal 
close packing. Grossman et al. employed free volume argu- 
ments in the vicinity of the close packing, and suggested an 
interpolation between the hexagonal-packing limit and the 
well-known low-density relations. The resulting CRs 11311 
read 



n c + n u.n(al + d) 2 T l l 2 
p = nT , K = 
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Here I is the mean free path, which is given by an interpola- 
tion formula 11311 : 
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where a = 1 — (3/8) 1 ' 2 . The CRs by Grossman et al. include 
three adjustable parameters of order unity: a, J and ji, where 



the latter drops out from the steady-state problem. Grossman 
et al. determined the optimum values a = 1 .15 and y = 2.26, 
from a comparison between MD simulations of a system of 
inelastic hard disks in a rectangular box, driven by a thermal 
wall, and numerical solutions of the hydrostatic equations 
(0 in rectangular geometry. We adopted the same values of 
a and y for the circular geometry. 

Employing Eqs. (|4]) and (0 we can reduce Eqs. (O to 
a single equation for the rescaled inverse density z(r,(j>) = 
n c /n(r, 0). In the rescaled coordinates r/R — > r the circle's 
radius is 1, and the governing equation becomes 
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is the hydrodynamic inelasticity parameter introduced in Ref. 
120]. As the total number of particles N is fixed, z _1 (r, 0) 
satisfies a normalization condition: 



/ d<j> dr rz~ l (r, </>) = nf, 
Jo Jo 



where 

2n \R 



(10) 



(11) 



is the average area fraction of the particles. 

Azymuthally-symmetric solutions of Eq. @, z(r,(p) — 
Z(r), are described by the following ordinary differential 
equation: 
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r dr 
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dZ 
dr 



= AQ{Z), 



while the normalization condition (1 1 Ob becomes 
l 

Z~ x (r)rdr = f/2. 



(12) 



(13) 



Equations ( 112b and (113b . together with the obvious boundary 
condition 



dZ 
dr 



0. 



(14) 



r=0 



form a complete set. The hydrostatic density profile is com- 
pletely determined by two scaled parameters A and /, and 
can be found numerically as we will show shortly. 
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Are the azimuthally symmetric states stable with respect 
to small perturbations? We performed marginal stability anal- 
ysis to find out whether there are steady-state solutions with 
broken azimuthal symmetry, z(r, </>) that bifurcate continu- 
ously from an azimuthally symmetric solution Z(r). Under 
additional assumption that the possible instability of the az- 
imuthally symmetric state is purely growing (that is, not os- 
cillatory), the marginal stability analysis yields the instabil- 
ity borders. The marginal stability analysis goes along the 
same lines as that developed for the rectangular geometry 
El 0; US HH. Let us search a steady-state close to an 
azymuthally-symmetric state: 



z(r,<p) =Z(r) +e <P k {r) sin(fc</>) , £<1. 



(15) 



As z{r, <j> + 2k) = z(r, </>), k must be an integer which can be 
chosen to be non-negative. Substituting Eq. dl5l) into Eq. (0I 
and linearizing around the azimuthally symmetric state Z(r), 
we obtain a linear eigenvalue problem, where k = k(f,A) 
plays the role of the eigenvalue: 



k 2 , AQ'(Z) 



F(Z) 



C* = o, 



(16) 



where 5fe( r ) = < ^'k{ r )F{Z). Equation dT6b is complemented 
by the boundary conditions 



6(0) =0 and &(1)=0 



(17) 



and can be solved numerically. Let us ignore for a moment 
the quantization of the eigenvalue k and, while looking for k, 
assume that it is a (positive) real number. In that case, a nu- 
merical solution yields, for a fixed A, a curve k = k{f), see 
Fig. |2] for two examples. At a fixed k, the azimuthally sym- 
metric state is unstable within an interval of area fractions. 
The foots of one such curve corresponds to the (hypotheti- 
cal) case when k tends to zero (that is, the azimuthal wave- 
length tends to infinity). The instability interval becomes 
narrower when k is increased, and it shrinks to a point at a 
maximum k max , signalling that a density modulation with a 
sufficiently short azimuthal wavelength should be stable for 
all /. For example, when A = 5 x 10 4 , the marginal stability 
curve k — k(f) has its maximum at k max » 0.38 (see Fig. [2j 
that is less than unity. Going back to the physical case, where 
k is a positive integer, we see that all the values for k, deter- 
mined from the eigenvalue problem, are unphysical. That is, 
there are no solutions to the eigenvalue problem that would 
satisfy the boundary conditions and the quantization condi- 
tion k = 1,2,3 We observed a similar behavior for values 

of A up to 10 8 . Though k max increases with A, the increase 
is extremely slow: slower than logarithmical (see the inset 
of Fig. [2]». These numerical results strongly indicate that the 
azimuthally symmetric states are stable with respect to small 
perturbations. This is in marked contrast with the presence of 
bifurcating states with broken symmetry in similar settings 
of a granular gas driven by a thermal wall, but in rectangular 
fl6l: fl7l: I20I: [l9l: l23h and annular [8] geometries. The absence 
of the bifurcating states with broken symmetry in the circu- 
lar geometry gives a natural explanation to the persistence 
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Fig. 2 The marginal stability curves k = k(f) for A = 5 x 10 4 and 
A = 10 5 , as indicated. The shaded area denotes in both cases the (un- 
physical) instability region that is obtained if one ignores the quantiza- 
tion of k: k= 1,2, The inset shows k max as a function of In (A ) . 



of circle-shaped cluster shapes as observed in our MD sim- 
ulations. 

Now let us compare the azimuthally symmetric density 
profiles, found from our hydrostatic calculations, with the 
results of MD simulations. As noted above, in the MD sim- 
ulations the coefficient of normal restitution £ = 0.888 was 
fixed, while the number of particles N was varied. Accord- 
ingly, the scaled parameter A = 9980 was fixed but / was 
varied. Figure [3] shows a comparison of the hydrostatic ra- 
dial density profiles with the radial density profiles obtained 
in MD simulations, for four different values of N. The MD 
profiles were averaged over 1000 snapshots. It can be seen 
from Fig. [3] that for N = 1716 the theory overestimates the 
density in the center of the box. This is expected as, for rel- 
atively small N, the maximum density is considerably less 
than the close-packing density, and the accuracy of the CRs 
by Grossman et al. is not as good. For larger N the agree- 
ment rapidly improves. We also computed the density pro- 
files using another set of semi-empiric CRs: those obtained 
in the spirit of Enskog theory ifbUl . One can clearly see from 
Fig. [3] that the Enskog-type CRs predict unphysically high 
densities in the cluster, and is also less accurate at interme- 
diate densities. Figure |4] compares, at different N, the max- 
imum radial densities predicted by hydrostatic theory with 
the CRs by Grossman et al. and those obtained in the MD 
simulations. The agreement is very good for the densities 
approaching the close packing density. 



4 Fluctuations 

Now we turn to stationary fluctuations of the radial coordi- 
nate r m (t) of the center of mass of the system. The radial 
PDF P(r m ,t) is normalized by the condition 



2n / P(r m ,t)r m dr m = 1 . 



(18) 
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Fig. 3 Scaled density n/n c versus scaled radius r/R as observed in 
MD (circles) and predicted by the hydrostatic theory with the CRs by 
Grossman et al. 11311 (solid lines) and with the Enskog-type CRs 1 14] 
(dashed lines) for four different values of the total number of parti- 
cles A'. For the rest of parameters please see the text. The dotted lines 
indicate n/n c = 1. 




MD simulations 
Hydrostatics 



_i_ 



_i_ 



1000 2000 3000 4000 5000 6000 
N 



0.08 



to 



0.06 




0.04 



0.02 



Fig. 4 The top panel depicts the maximum scaled density n max /n c as 
a function of N as predicted by the hydrostatic theory with the CRs 
by Grossman et al. 11311 and observed in MD simulations. The bottom 
panel shows the scaled standard deviation of the center of mass posi- 
tion as a function of N, obtained in MD simulations. 
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Fig. 5 The logarithm of the radial PDF P(r m ) of the center of mass 
position versus (r„,/R) 2 , for four values of N. 



where we have returned to the dimensional coordinate. Typi- 
cal MD results are presented in Fig. [5] which shows logP(r m ) 
versus (r,„/R) 2 , for four values of ,/vQThe observed straight 
lines clearly indicate a Gaussian. This finding suggests a 
Langevin description of the macro-particle. Consider a macro- 
particle, performing an over-damped motion in a confining 
harmonic potential U (r m ) = kr^/2 and driven by a (delta- 
correlated) discrete-particle noise r](t). The Langevin equa- 
tion for this problem reads 



hf m + kr m = T](t), 



(19) 



where h is the damping rate, (rf(t)ri(t')) = 2/zF<5(f-f'), and 
F is an effective magnitude of the discrete particle noise. In 
the limit F — > 0, the (deterministic) steady-state solution of 
Eq. d 19b is r m = f m = 0: the macro-particle at rest, located at 
the center of the box. At F > 0, the steady state PDF is (see, 
e.g. Ref. JH) 



P(r m ) 



1 



exp 



(20) 



where a 2 = T /k, and the normalization constant is com- 
puted for a <C R. The variance a 2 is the ratio of F (a charac- 
teristic of the discrete noise) and k (a macroscopic quantity). 
Unfortunately, the present state of theory does not enable us 
to calculate either of these two quantities. An important in- 
sight, however, can be achieved from the A'-dependence of 
a, obtained by MD. In analogy with equilibrium systems, 
one might expect the relative magnitude of fluctuations to 
decrease with increasing N. Surprisingly, this is not what 
we observed, see the bottom panel of Fig. [4] One can see 

1 For larger values of N the density in the center of the container 
approaches the close packing (see Fig.O, and the particle collision rate 
in the high-density region almost diverges. As a result, the event-driven 
simulation progresses extremely slowly, and the simulated real time is 
not large enough to allow for good statistics needed to determine a 
reliable PDF. 
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that o(N) approaches a plateau, that is fluctuations persist at 
large 



5 Discussion 

We employed a simple model system [9] to investigate the 
density profiles and fluctuation properties of dense clusters 
emerging in driven granular gases. The density profiles, ob- 
tained in our MD simulations, are well described by hy- 
drostatic equations which employ the constitutive relations 
suggested by Grossman et al. 111311 . The good performance 
of these constitutive relations is in agreement with previ- 
ous results in rectangular geometry, with and without gravity 
fl3t l22h . Marginal stability analysis yields a natural expla- 
nation to the circular cluster shape, observed in the MD sim- 
ulations. Finally, the observed Gaussian fluctuations of the 
center of mass suggest an effective Langevin description in 
terms of a macro-particle in a confining potential of hydro- 
dynamic nature, driven by discrete particle noise. The fluc- 
tuations persist as the number of particles in the system is 
increased, and this surprising finding awaits a proper theo- 
retical interpretation. 
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